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A CLASS OF MONOTONIC SHOCK-CAPTURING DIFFERENCE SCHEMES 
0. Introduction 

A broad class of gas-dynamic problems can be treated within the 

framework of ideal gas theory, i.e., assuming the gas is inviscid and 

thermally nonconductive. The differential equations of gas dynamics 

follow from the laws of conservation of mass, momentum, and energy 

in integral form assuming continuous differentiability of the fluid 

variables, and constitute a system of quasi-linear hyperbolic equations* 

It is well known that in general these may posess discontinuous solutions 

even when smooth initial data are prescribed^. Physically, the presence 

of discontinuities in the solution typically signals the appearance of 

a shock wave. There are two possibilities for correctly describing 

discontinuous flows within the framework of ideal gas theory. The 

first consists of breaking up all the regions In the problem into 

subregions of smooth flow, which are described by the differential 

equations of gas dynamics, while the discontinuities (boundaries of the 

subregions), are described by conservation conditions. It is important 

to distinguish weak discontinuities (discontinuities in derivatives of 
the fluid variables), tangential discontinuties, and finally, strong 

discontinuities (shock waves). The second approach consists in utilizing 

the conservation conditions in integral form, which allows discontinuous 
solutions. For historical reasons (the differential equations were 

derived earlier), these are referred to as "generalized," iu contrast 
to the continuous classical equations. We note that consideration of 

generalized solutions extends the class of possible solutions, and it 
is therefore necessary to make use of additional considerations in order 

to determine the correct solution. For example the requirement that 
entropy not decrease at a discontinuity permits us to exclude 
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rarefaction shocks in the flow of a perfect gas (from a formal point 

of view such a discontinuity is unstable). 

In numerical modelling of gas motions in the presence of shock waves 

both of the above approaches are possible. In the former, the physical 

jump conditions on the discontinuities dividing the regions of continuous 

flow are used both in order to obtain the conditions connecting the fluid 

variables on the two sides of the discontinuity and to determine the 

motion of the computational mesh points following the discontinuity. 

Changes in the topology of the lines of discontinuity can be followed 

either using a moving curvilinear mesh or by means of an algorithm based 

on a fixed mesh. The first method is most completely worked out in the 

2 

paper by S. K. Godunov et al. while the second has been developed by 

. , 3,4 

Moretti et al. 

The shortcoming of both approaches is the inhomogeneity of the 
resulting difference schemes and consequent complicated structure of tb> j 
numerical algorithm. Both approaches, as a rule, demand an a p riori know¬ 
ledge of the flow pattern during the calculation. 

When it is impossible to know in advance the flow pattern t 6r when 
it is changing qualitatively with time, it is more convenient to make 
of "shock capturing” difference schemes, amounting to difference approxi¬ 
mations to the integral form of the conservation laws for every compute-* 
tional cell. The form of the difference equations does not depend on 

the character of the flow or the position of the possible shocks, and 

1 5 

therefore such schemes are homogeneous. * [Apparently homogeneous 11 

jC 

means universal, i.e., not problem-dependent.—DLB) In this approach 

* * 

shocks appear as regions of abrupt variation of the fluid quantities. 
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Shock-capturing schemes can also be based on the differential equations. 
For this purpose small corrections are introduced in the equations, 
typically of nonlinear form, analogous to physical viscosity. The 
equations assume parabolic form and permit a smooth solution which 
tends to the solution of the original system as the tf artificial ,f vis¬ 
cosity vanishes. However, difference analogs of the conservation laws 
are satisfied, in this technique, generally speaking, only approximately. 
In difference schemes based on the integral form, the mass, momentum, 
and energy conservation laws are sati^-'ed for every computational cell 
to roundoff. Such schemes are describable in terms of fluxes across 
boundaries of the computational cells, and so the conservation laws 
for each computational region are algebraic consequences of the con¬ 
servation laws for cells, i.e., the schemes are conservative,^ 

We note that the first approach, as a rule, distinguishes only 
simple shocks, while the other types of discontinuity are calculated 
as in the shock-capturing method. A difference scheme employed for 
solving practical problems is required to be accurate, that is, it 
must be able to describe flows on comparatively coarse meashes; it must 
be economical; and it should be simple to employ. The accuracy of the 
scheme for calculating smooth flows is determined by the order of the 
approximation. In the neighborhood of a shock, the change In the fluid 
variables is comparable with their magnitude, so that the concept of 
the order of the approximation becomes meaningless. It is well known 
that schemes of the first type lead to smearing of the shock because 
of the strong numerical viscosity. Schemes of second and higher order 
give rise to significantly less smearing of shock runs, but are typically 









nonmonotonic. One of the principle reasons for the appearance of un- 
physical oscillations in the neighborhood of a shock is the nonvanishing 
dispersion of the difference scheme, which leads to misrepresentations 
of the form and speed of propagation of physical disturbances, especially 
at short wavelengths. Nonmonotonicity may also be produced by the non¬ 
linearity inherent in real problems.® 

Artificial viscosity is introduced to suppress unphysical oscil¬ 
lations. Most methods for achieving this, however, do not yield a 
monotonic scheme, and consequently make the localization of shocks worse 

and strongly smear large gradients* This problem becomes especially 

* 

severe in calculations involving strongly shocked flow arising from 

nonsteady interaction of shock waves with one another and with obstacles* 

Recently a number of methods using nonlinear limiters of various sorts 

7-14 

to filter unphysical oscillations have appeared. Many of these are 

highly effective, so that the assertion that "good schemes of first and 

2 

of higher order smear shock fronts in practically the same degree” 
cannot be accepted as correct. 

It is convenient to choose an actual difference scheme in several 
stages. First the dispersion and dissipation properties are studied 
in the linear approximation using the method of differential approxi¬ 
mations or harmonic analysis for the simplest transport equations. It 
is essential to test the schemes selected through linear analysis on 
simple one-dimensional problems having an exact solution. The final 

stage of development of any difference scheme must be its verification 

< 

on model two-dimensional problems. 





In the present paper a nonlinear conservative smoothing procedure 

is proposed in order to achieve monotonicity In explicit schemes of 

higher order.' Although this smoothing procedure is applicable to 

arbitrary schemes, the underlying transport algorithm is taken here to 

15 

be second-order MacCormack differencing. The results of the linear 
analysis of the dispersion and dissipation properties of this and a 
number of other widely used schemes are presented in Ref. 16. Below, 
the second stage in the choice of an optimal scheme is considered in 
detail. This is conveniently divided into two parts. The first is an 
investigation of the solution of the linear advectlon equation, which 
facilitates an understanding of the nonlinear properties of the scheme. 

The second part is a calculation of the flow resulting from the evolution 
of a one-dimensional gas dynamic shock. In conclusion, as an illus¬ 
tration of the possibilities of the proposed numerical method, a 
three-dimensional problem involving the interaction of a plane shock 
wave with an obstacle is discussed. 
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1. Construction of second order monotonic difference schemes 


An effective device for constructing monotonlc difference schemes 

having second-order accuracy for smooth flows was presented by Boris 
11 12 

and Book. * It consists of two stages. In the first a large numeri 
cal diffusion is Introduced into the solution, which guarantees the 
monotonicity of the scheme. In the second stage this diffusion is 
canceled wherever doing so does not introduce new unphysical extrema or 

accentuate existing ones. By means of this approach a family 
of FCT (Flux-Corrected Transport) methods was contetructed. 

We consider the simple example of the Cauchy problem for the 

* 

linear advection equation 


where 



v 


3f 

3x 


0 


( 1 . 1 ) 


f(0,x) *= f Q (x). 

The exact solution has the form f(t,x) = fg(x ~ vt). Let f^ be the 
value of the function at the i* . grid point of a uniform mesh on the nth 
time level; let i 4* T be the operator which propagates the solution 
from the nth to the (n+l)st level. The form of T depends on the 
particular difference scheme used. We define a diffusion operator D by 


Df i " Vl/2 ' Vl/2 " Q6f i+l/2 ~ Q6f i-l/2 

where 

6f i+l /2 * f i+l “ f r 


( 1 . 2 ) 






Then the first (diffusion) stage can be represented in the fora 

* (1+T+D)f°. 

Now we cancel the diffusion in such a way that in the resulting solu- 
tion no new extrema In comparison with appear, and those already 
present are not amplified. For this purpose in the antidiffusion 
operator A, given by 


M i (<? i+l/2 “ ^i-1/2* * 


we limit the fluxes = ^*±+ 1/2 accordin 8 t0 the formula 


*1+1/2 " S * Bax t0 * 6? i-l/2‘ I'Pi+l/zl 


S ^i+3/2 1 *’ 


(1.3) 


where S = sign cp.+l/2* 

The transition from the nth to the (n+l)st time level has the form 


f” +1 = (1+A)7^ = (1+A)(l+T+D)f“ . (1.4) 


This technique—the explicit cancellation of the diffusion—leads to 
retention of some of the diffusion in smooth regions [i.e., where the 
limiter (1.3) does not operate and A - D], even when T 5 0. Conse- 

4 * 

quently Boris and Book proposed two other algorithms: 
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... a. 
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,, phoenical ,, 

fj +1 - [(1+A)(1+T) + D]fJ (1.5) 

and 

implicit 

f” +1 = (l+Aj^d+T+D)^, (1.6) 

which permit complete cancellation of the diffusion when T = 0 in 

e 

smooth regions. We note that using the algorithm (1.4)-(1.6) leads to 
substantial improvement in the results in comparison with widely used 
schemes. 

In the SHASTA^ scheme the coefficient Q *» 1/8. Subsequently 

Boris and Book investigated in detail the influence of the value of 

magnitude of Q on the dissipative properties of this scheme and deter- 

12 

mined the optimum dependence of Q on the Courant number. However, 
practically speaking, this optimization is meaningful only for a linear 
equation with constant coefficients. For the equations of gas-dynamics 
as calculations reveal, good results are obtained for Q between 
1/10 and 1/6. 

One can point out three factors responsible for the success of the 
method of Boris and Book. First, the diffusion stage, as shown by 
linear analysis, substantially improves the dispersion properties of 
the scheme. Second, the diffusion is introduced in a conservative 
manner, and, finally, it is done so as to permit excellent localization 
in the neighborhood of the shocks. Evidently it is precisely these 
last two properties which are most important. In Ref. 13 it was . 



proposed to interpret the method of Boris and Book as a nonlinear con¬ 
servative smoothing procedure and to introduce diffusion at the (n+l)st 
time level according to 


f* - (1+A)(l+D)(l+T)fj (1.7) 

or 

f" +1 - (1+A+D)(l+T)f“ <1.8) 


This procedure fails to improve the dispersive properties and falls 
somewhat short in quality of the results of the algorithms (1.4)-(1.6). 
However, in contrast to the original method of Boris and Book, which is 
essentially one-dimensional in character and which in two- 
dimensional calculations can only be applied via the method of coordi¬ 
nate timestep-splitting of the the difference equations, smoothing in 
the fora of Eqs. (1.7)-(1.8) is easily incorporated in essentially 
two-dimensional schemes. For this purpose, after each timestep successive 
smoothing operations are carried out with respect to each of the coordi¬ 
nates. Such smoothing has turned out to be effective in solving a 
whole series of problems. 

In the present paper a simpler smoothing procedure is proposed, 
having a local character while retaining the conservative property of 
the difference scheme. A constant diffusion is introduced only in 
regions which are not monotonic, i.e., where the solution has a numerical 
ripple. As a result of this, just as in the method of Boris and Book, 
there is a certain amount of smearing of physical extrema, but for.the 

4 ‘ 

most part only oscillations produced by the numerical scheme are 

9 









removed. The conservative property is assured by introducing the 
diffusion in the form of fluxes across the cell boundaries. Using 
the notation introduced above, we can write the smoothing at the nth 
time level as 


f" +1 - (1+T)fj + D*fJ, (1.9) 

and on the (n+l)st time level in the form 

t 

( 1 . 10 ) 

D f i ■ %H./2 ~ *1-1 fV 

*1+1/2’ 6f 1+1/2 6f 1+3/2 < 0 or 1+1/2 6f 1-1/2 < 0 

(l.H) 

I 

0 otherwise. 

In order to code this prescription it is convenient to search the entire 

array f., recording the boundaries across which the flux is nonvanishing, 
x 

and then in a second pass calculate the nonzero fluxes. 

This procedure for introducing artificial viscosity runs substan¬ 
tially faster than the method of Boris and Book, while yielding com¬ 
parable results. Note that, in contrast with the algorithms (1.4)-(1.7), 
smoothing according to the prescription (1?8)-(1.10) in regions where 
the fluid quantities vary monotonically in general leaves the solution 
unchanged. 


where 


*1+1/2 


f^ 1 - (1+D*)(l+T)f* . 

In these expressions 
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2. Solution of the linear transport equation 

In this and the following section the proposed numerical method 

is compared with a number of different schemes, using a one dimensional 

22 

model problem having an exact solution. The schemes of S. T. Godunov , 

15 23 

HacCormack , V. V. Eremin and Yu. M. Lipnitskii , of first, second 

and third order, respectively, are considered and also the non-linear 

12 

scheme of Boris and Book, SHASTX , embodying the algorithm (1.5). These 
schemes are widely employed and have performed well in practice; therefore 
the present analysis permits one to carry out an objective evaluation of 
the proposed method. 

Let us approximate equation (1.1) on a uniform mesh with Ax * 1. 

The computational mesh contains a total of 100 cells. On the boundaries 

The initial conditions are given in 

x > L 
1 £ x £ L 


periodicity conditions are imposed 


the form 


f(o 


(0.5, 

» x) " (cp(x), 


where L « 20 and 


cp(x) = 2 


( 2 . 1 ) 


or 

<p(x) - 0.5 4* 0.075x. (2.2) 

In Fig., 1 are shown the results of the calculations after 800 

time steps with Courant number v« 0*2, i.e. for a total running 

time t - Vt/L, after which the initial square wave (2.1) of width 20 . 

cells has been carried across 160 cells. The curves (al), (a2), (a3), 

correspond to the schemes of Godunov (which coincides in the linear 

case with one-sided differencing), MacCormack (equivalent to the scheme 

«* 

of Lax and Wendroff), and Eremin-Lipnitskii; (b), to the SHASTX scheme; (c). 
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(d), (e) correspond to the monotonic MacCormack scheme with smoothing 
according to algorithms (1.7), (1.10), and (1.9). The best results 
were obtained with SHASTX. The SHASTA** algorithm and the monotonic 
MacCormack scheme with the smoothing (1.9) were very slightly inferior 
to this. The smoothing procedures (1.7) and (1.10) give worse results, 
which are comparable, however, with the results obtained from the third- 
order scheme. The first-order and second-order schemes yielded 
unsatisfactory approximations to the solution. 

Since the solutions of the linear advection ( equation using the 

« 

algorithms (1.7) and (1.10) were inferior, however slightly, to the algorithm 
(1.9), in the remainder of this section we will present the results of 
calculations carried out using the latter. We must point out, however, 
that these other algorithms will be analyzed below, inasmuch as they 
possess definite advantages in solving real gas-dynamic problems. The 
results of the present section are applicable in problems where linear 
equations are solved, for example, in certain meteorological problems. 

Figure 2 illustrates the ”flexibility" of the various schemes in the 
sense of Ref. 2, i.e., the capability of solving problems in which the 
Courant number varies greatly over the computational region. Here 
profiles are shown at t * 3, obtained for Courant number v* 0.6 
(curve 1) and v= 0.2 (curve 2) using first-order (a), second-order (b), 
third-order (c) and the monotonic MacCormack (d) schemes * The 
superiority of the latter is indubitable. The calculation was not 
carried out with SHASTX, which requires that condition v< 0.5 be 
satisfied**. « 

m 

As we remarked in Section 1, the smoothing which removes numerical 
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oscillations smears out physical extrema. Consequently, the test 
problem using initial conditions in the form of a triangular profile 
(2.2) constitutes a stringent test of a scheme. Figure 3 shows 
profiles at time t ** 1 (curves 1) and t * 20 (curves 2) for the same 
schemes as in Fig. 2, if v- 0.2. The first-order scheme practically 
"forgets" the Initial conditions; the second-order scheme qualitatively 
changes the shape of the pulse as time elapses, while the thld-order 
scheme preserves the "height" of the distrubance well. The monotonic 
MacCormack scheme is only negligibly inferior to SHASTX (the dashed lines 
in Fig. 3d). 


I 





3* The evolution of a gas-dynamic shock 

The problem of the evolution of a gas dynamic shock Involves all the 
characteristic features of one-dimetu'ional ideal gas flow: creation and 
propagation of a shock wave, contact discontinuity and rarefaction fan. 

By comparing the results of the calculation with the analytic solution, 
ve can check both the accuracy with which the characteristic flow regions 
are transported and the speed with which they develop-, A comparison 
of the various schemes mentioned above will be presented, based on this 


problem. 

As was stated earlier, unphysicaloscillations in the neighborhood of 

shocks and steep gradients appear because of. the nonvanishing dispersion 

associated with a difference scheme, the nonlinearity of the gas-dynamic 

24 

equations, and also for several other reasons • Besides this, it is 

possible that the various nonlinear smoothing operators may Interact 

with one another nonselfconsistently. 

The system of equations has the following form: 


f p| dx + f pul' 
•'x, ’t. J t, *x. 


dt = 0, 


*2 t 2 t 2 | x 2 

I pu J dx + f (pu 2 + p)| 
*' x 1 *1 **1 %1 


dt » 0. 


C 2 \ 2 C 1 r 

J e I dx + J (e + p )u I 

*v t, X, 


dt * 0. 


Here e *= p/(y-l) + pu /2, where y is the adiabatic index. We consider 

< 

the evolution of a shock with initial pressure ratio p^/p^ “ 2 and 
density ratio pj/p^ “ 1» and y * 1,4. In Fig. 4 are shown density 
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m 


profiles obtained using the Godunov (a) , Eremin-Lipnitskii (b), 

MacCormack (c) and SHASTX (d) schemes. The exact solution Is 
indicated by the heavy continuous line and the numerical results 
up to t B 25 by a fine continuous line. The first-order scheme, as 
one might expect, badly smears the shock front and contact discontinuity 
and represents the rarefaction fan poorly. The Eremin-Lipnitskii and 
MacCormack schemes are nonmonotonic and yield oscillations in the 
neighborhood of the shocks. SHASTX, which represents the shock front 
very well, gives rise to small undershoots in the density in the neighborhood 
of the contact discontinuity and rarefaction fan. For this case the 
SHASTA scheme, as described by Boris and Book, gives practically the 
same results as does SHASTX. All of the schemes accurately represent 
the pressure profile between the shock wave and the rarefaction fan. 

In Fig. 5 are shown results of calculations using the monotonic 
MacCormack scheme with various forms of smoothing. Fig. 5a corresponds 
to algorithm (1.7). One shortcoming of this approach lies in the presence 
of undershoots in density in the neighborhood of the shocks and rarefaction 
fans. Very similar results are obtained when one applies the smoothing 
(1.9) to all of the equations in system (3.1) (Fig. 5b). 

In the present work it is proposed to carry out self-consistent 
smoothing of the density, momentum, and energy in regions where 
nonmonotonicity in the density profile exists, using the procedures 
(1.9) and (1.10) described above. The justification for this lies in 
the fact that every stable one-dimensional shock has a density 
discontinuity. In particular, this self-^consistent smoothing 
introduces no changes in the velocity at a contact discontinuity, since 
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the density and iaosenturn are smoothed in a "related” manner* It also 
somwehat reduces the expenditure of machine time, since In order to 
determine nonmonetonicity, one array of numbers instead of three is 
searched [Cf.(l.ll)}. We note that this type of smoothing cannot be 
implimented using algorithms (1*4) - (1*8) because of their two-stage 
character• 

The density profiles obtained using self-consistent smoothing at 

th ' 

the n and (n+l) st levels are shown in Figs. 5c and 5d. These 

results are clearly superior to all those shown previously. It is 

interesting to observe that the difference in the results obtained by 

t 

smoothing according to the algorithms (1*9) and 1.10) is much smaller 
than in the case of linear advection (Fig* 1) . This is indicative of 
the decisive role played by nonlinear effects in numerical solutions 
of the gas-dynamic equations. 

We pause to draw attention to some peculiarities in the solutions 
we have found. Most of the calculations were run with At « 0.5. 

SHASTA, SHASTX, and the MacCormack scheme with self-consistent 
smoothing, however, yielded a highly oscillatory solution. The 
oscillations are a consequence of the nonselfconsistency of the smoothing 
and the interaction between the resulting errors in the various fluid 
quantities. A separate paper will be. devoted to a detailed analysis 
of the nonlinear properties of the different schemes, the nonmonotonic5 
resulting from various causes, and ways of contending with the unphysical 
oscillations that result. We therefore will not go inot the reasons for 
the appearance of computational ripples in detail, but content ourselves 
with remarking that they vanish as the time step is reduced. The 1 
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solutions shewn in Fig. 4c and 5b were obtained with a time step 
t - 0.25. 

We now compare the respective running times. Taking as the unit of 
time that required for a single step according to the MacCormack scheme, 
the corresponding running times for the other schemes ate as follows: 
monotonic MacCormack with self consistent smoothing, 1.4; with 
nonselfconsistent smoothing according to algorithm (1.9), 1.5; using 
algorithm (1*7), 2; for the Godunov scheme 2 - 2.5, depending on how 
the "large" quantities are calculated; for SHASTA and SHA;;TX, 3; and 
for the Eremin-Lipnitskii scheme, 2.5. 

We next consider an example in which a strong shock wave is formed. 

The initial parameters are taken to be the same as in the paper of Boris 
11 

and Book (P2^P^ m 480, p 2^ p l " Y “ 5/3, At: « 0.02). The Eremin- 
Lipnitskii and MacCormack schemes in the absence of ai> additional 
artificial viscosity do not permit calculation of such a Oov, SDAvrA 
and the monotonic MacCormack scheme with smoothing according to 
algorithm (1.7) give approximately the s^ame results as SHaSTX, We 
• therefore make further comparisons of the results obtained, using the 
method of Godunov, SKASTX, and mcnotonir. MacCormack scheme with 
consistent smoothing given by (1.9) , In Figs. 6 and 7 a) c shown density 
and pressure profiles obtained using these nebr-men at tim* 1 ^ 2 (deshrd 
lines and t *» 4 (fine continuous lines), A vertical dashed- dotted lint 
is used to show the initial position of the shock. The vjvt uec: and 
shortcomings of the numerical methods under investigat ion nre 
clearly shown i** Table 1, in which the values of density and relative 
errors at time t » 4 are shown. In the left column of the table is shown 
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the number of computational cells. The Godunov scheme strongly 
smears the shock front and contact discontinuity and approximates the 
exact solution poorly in the region of the rarefaction fan. SHASTX 
is someiv’hat better than the monotonic MacCormack scheme in representing 
the shock wave and rarefaction fan, but somewhat inferior to It in 
representing the contact discontinuity* In addition, SHASTX gives 
incorrect values In the neithborhood of the Initial shock* 

The mono tonic MacCormack scheme and SHASTX have good dynamical 
properties. They permit relatively fast determination of the true flow 

pattern (Cf. the Godunov scheme). This last property is especially Important 

* 

in solving nonstationary problems. 
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4. Calculation of three dimensional interaction of a shock wave with 
an obstacle* 

The possibilities of the proposed method are illustrated below for 

the example of reflection of a plane shock wave from a obstacle of 

finite width. The initial stage of the interaction precess is 

investigated, in which the reflected shock wave begins to form and the 

pressure loading on the object is a maximum. We remark that this problem 

is of interest in its own right, inasmuch as experimental determination of 

nonstationary three-dimensional flows encounters serious difficulties. 

The problem is solved in Cartesian coordinates. The mesh is 

uniform. Timestep splitting was employed in solving the difference 
25 26 

equations 9 which permitted a substantial saving in running time. 

The conservation laws in Eulerian coordinates have the form 


J" p J ^ dv + J 2 (j) pV*n 

v C 1 t l s 

/pv | 2 dv+ / 2 J) (pW + p) 

V s 


f | fc 2 f t 2 r - - 

§ e I dv + I (D (e + p) V • n 
_ _ 1 t _ s 


ds dt * 0, 


n ds dt 8 0, 


ds dt 13 0 


We denote by (At) the operator which performs a step in the direction a: 


f 


n + 1/3 
ijk 



(At) f 


a 

ijk. 




L is defined in the following manner: 
a 
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T ljk ' f Jk - i (f A,J + °’ Kks- r ijk), 


n +i/3 m 1 r f n .V At ~ 

f ijk 2 Lf ijk + f ijk " Ah CF ijk " F i-q,j-r,k-s' J » 


where 


/p 1 



f° 

pu J 

/puj 


[q 

P U y 

f “K 

u + 
a 

X 

[ou! 

\ PU 2 


[ S 

In 

l e 1 

H 

[ VI 


P» 


q * S , r » 6 , 

n xa ya 


za. 


The smoothing is carried out according to the algorithm (1.9) • The 
transverse momentum components were smoothed independently. 

The complete transport operator will look as follows: 

L(2At) - L (At) L (At) L (At) L (At) L (At), 
y x z x y 

Written thus, even if the split operators L o (At) do not commute, the 

27 

difference scheme has second-order accuracy over a double tlmestep 
The stability condition for the split equations is 


At £ min (At ), a * x,y,z, 
a a 


where 


At £ min (~rr 
a ijk “a 48 




and a is the speed of sound. This condition turns put fco he 
considerably less stringent than for the unsplit operators: 


At £ min ( —— 

I 


Ah 


|v| + a/3 
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Since At ( x is the direction of propagation of the shock wave) 

is usually much smaller than At and At , it makes sense to carry out 

y z 

the calculation in the y and z directions with a double timestep. Thus 
using 


L Q (2At) - L a (At) L a (At), a = y,z. 


we finally obtain 


n+2mAt 


ijk 


L (At)[L(At) L (2At) L (At) L (2At)] 

Jr A £* & Y 


m-1 


n 


. L (At)L (2At)L (At)L (At)f . 
x z x y xj K 


Use of this algorithm instead of 


n+2mAt n 

f « L n (2At)f 
ijk 


permits reduction of the computational volume by a factor of 1.5. 

In the test the computationalmesh consisted of 20 X 20 X 15 cells. 

Calculation of a single case up to t = 1.2 - 1.5 required about 2 hours of 

time on the M4030. Here t « tV /h. where V is the speed of the 

amp amp r 

reflected shock wave which arose as a result of the interaction of the 
incident wave with the square object,, and h is the step height. 

A typical flow pattern is shown in Fig. 8. At t * 1.1 the isobars 
are shown in (a) with separation Ap = 10.0,and the lines of constant 
Mach number (b) at intervals AM - 0.2. The Mach number of the incident 
shock wave M =5.0 and the ratio of width to height of the step was 
t/h ® 5.2. The adiabatic index was y « 1.4 and tha pressure in the 




w 


iiL. j. min 




undisturbed gas was P^ = 1-0. It is clear that the method permits 
excellent resolution of the shock wave and the flow in the neighborhood 
of the corner points. 

In Pig. 9 are shown the isobars in two cross sections x « const* 
located equidistantly and parallel to the boundary of the step at a 
distance 0.3 h. The ratio i/h * 2.0 in (a), 3.6 in (b), and 5.2 in (c). 

In Fig. 9a are shown the isobars for the plane case (i/h * «) • 

Figure 10 illustrates the influence of the ratio i/h on the bending 
of the shock wave in the xy plane. Curves 1,2,3 ^correspond respectively 
to l/h « 2/0, 3.6, and 5.2 (the bending is related to the bending for the 
plane case)• 

The authors thank Yu. P. Golovachev for his unstinting attention to 
this work and useful discussions, and also V.M, Goloviznin, G*L. Stenchikov 
and A.P. Favorskii for valuable suggestions. 
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Figure 4 





















a (Godunov method) 



c (monotonic 
MacCormack scheme) 
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U.S. Army Mobility Equip R&D CMD 
Fort Belvoir, MD 22060 

(CNWDI to Army Mat Dev 
& Readiness Command) 

Olcy Attn DRDME-WC (Tech Lib) 

Commander 

U.S. Army Nuclear & Chemical Agency 
7500 Backlick Road 
Building 2073 
Springfield, VA 22150 

(desires only Icy to Library) 
Olcy Attn Library 
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Commander 

David Taylor Naval Ship 
R&D CTR 

Bethesda, MD 20084 
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Mrs. M. Birkhead 
Code 5815.61) 

Clcy Attn Code L42-3 
(Library) 

Of ficer-in-Charge* 
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Civil Engineering Laboratory 
Port Hueneme, CA 93041 

Olcy Attn Code L51 S Takahashi 
Olcy Attn Code L51 R Odello 
Olcy Attn Code L08A (Library) 
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Naval Electronic Systems Command 
Washington, DC 20360 

Olcy Attn PME 117-21 
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Naval Facilities Engineering Command 
Washington, DC 20390 

Olcy Attn Code 09M22C (Tech Lib) 
Olcy Attn Code 03T 
Olcy Attn Code 04B 

Headquarters 

Naval Material Command 

Washington, DC 20360 

Olcy Attn MAT 08T-22 

Commander 

Naval Ocean Systems Center 
San Diego, CA 92152 

Olcy Attn Code 4471 (Tech Lib) 
Olcy Attn Code 013 E Cooper 

Superintendent 

Naval Postgraduate School 

Monterey, CA 93940 

(desires no CNWDI Documents) 
Olcy Attn Code 0142 Library 

Commanding Officer 

Naval Research Laboratory 

Washington, DC 20375 

(RD & RD/N Attr Code 1221 for 
& FRD Attn Code 2628 for) 

Olcy Attn Code 2627 (Tech Lib) 
Olcy Attn Code Boris/Book 


Commander 

Naval Sea Systems Command 
Washington, DC 20362 

Olcy Attn Sea-09G53 (Lib) 

Olcy Attn Sea-0351 

Officer In Charge 
Naval Surface Weapons Center 
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Silver Spring, MD 20910 
Olcy Attn Code F31 
Olcy - Attn H. Giaz 
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Naval Surface Weapons Center 
Dahlgren, VA 22448 

Olcy Attn Tech Library & 

Infor Services Brnch 

President 

Naval War College 

Newport, RI 02840 

Olcy Attn Code E-ll 
(Tech Service) 

Commander 

Naval Weapons Center 
China Lake, CA 93555 

Olcy Accn Code 2 66 C Austin 
Olcy Attn Code 3201 P Cordle 
Olcy Attn Code 233 (Tech Lib) 

Commanding Officer 
Naval Weapons Evaluation Facility 
Kirtland Air Force Base 
Albuquerque, NM 87117 

Olcy Attn R Hughes 

Olcy Attn Code 10 (Tech Lib) 

Office of Naval Research 
Arlington, VA 22217 

Olcy Attn Code 474 N Perron 
Olcy Attn Code 715 (Tech Lib) 
(Unci only) 

Office of the Chief of Naval 
Operations 

Washington, DC 20350 
Olcy Attn OP 981 
Olcy Attn OP 03EG 

Director 

Strategic Systems Project Office 
Department of the Navy 
Washington, DC 20376 

Clcy Attn NSP-43 (Tech Lib) 
Olcy Attn NSP-272 
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documents) 
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Olcy 
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Air Force 
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Attn NTE M Plamondon 
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Attn -R-Guice 
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Attn DEX 

Assistant. 

Chief of Staff 


Intelligence 

Department of the Air Force 
Washington, DC 20330 
Olcv Attn IN 

hi*' lid ic Missie Office/DE 
Air Force Systems Command 
Norton AFE, CA 92409 
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Air Force Systems Command 
Norton AFB, CA 92409 

(Minuteman) MNNX 
Olcv Attn MNNXH D Gage 

Deputy Chief of Staff 
Research, Development, & ACC 
Department of the Air Force 
Washington, DC 20330 
Olcv Attn AFRDQSM 

Deputy Chief of Staff 
Logistics 6 Engineering 
Department of the Air Force 
Washington, DC 20330 
Olcy Attn LEEE 
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Foreign Technology Division* AFSC 
Wright-Patterson AFB, OH 45433 
Olcy Attn NILS Library 

Commander 

Rome Air Development Center, AFSC 
Griffiss AFB, ICY 13441 
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Olcy Attn TSLD 

Stragetic Air Command/XPFS 
Department of the Air Force 
Offutt AFB, NB 68113 

Olcy Attn NRI-STINF0 Library 
Olcy Attn XPFS 

Department of Energy 
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P O Box 5400 
Albuquerque, NM 87115 
Olcy* Attn CTID 

Department of Energy 
Washington, DC 20545 

Olcy Attn OMA/RD&T 

Department of Energy 
Nevada Operations Office 
P 0 Box 14100 
Las Vegas, NV 89114 

Olcy Attn Mail & Records 
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Department of the Interior 
Bureau of Mines 
Bldg. 20, Denver Federal Ctr 
Denver, CO 80225 

Olcy Attn Tech Lib (Unci only) 

Director 

Federal Emergency Management Agency 
1721 I Street, NW 
Washington, DC 20472 

Olcy Attn Hazard Eval & Vul 
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Artec Associates, Inc. 
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Olcy Attn Library A830 

BDM Corp. 
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Olcy A:in T Neighbors 
Olcy At:n Cor r orate Library 

BDM. Corp. 

P 0 Box 9274 
Albuquerque, XM 87119 
Olcy Attn R Hensley 

Boeing Co. 
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Seattle, WA 98124 

Olcy Attn M/S 42/37 R Carlson 
Olcy Attn Aerospace Library 

California Research 6 Technology, Inc. 
6269 Variel Avenue 
Woodland Hills, CA 91364 
Olcy Attn Library 
Olcy Attn K Kreyenhagen 

California Research & Tech, Inc. 

4049 First Street 
Livermore, CA 94550 
Olcy Attn l> Orphal 


Calspan Corp. 

P O Box 400 
Buffalo, NT 14225 
Olcy Attn Library 

Denver, University of 
Colorado Seminary 
Denver Research Institute 
P 0 Box 10127 
Denver, CO 80210 
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J Wisotski 

EG&G Washington Analytical 
Services Center, Inc. 

P 0 Box 10218 
Albuquerque, NM 87114 
Olcy Attn Library 

Eric H. Wang 
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University of New Mexico 
University Station 
P 0 Box 25 

Albuquerque, NM 87131 
Olcy Attn N Baum 

Gard, Inc. 

7449 N Natchez Avenue 
Niles, IL 60648 

Olcy Attn G Neidhardt 
(Unci only) 

General Electric, Co 
Space Division 
Valley Forge Space Center 
P 0 Box 8555 
Philadelphia, PA 19101 
Olcy Attn M Bortner 

General Electric Co.-Tempo 
816 State Street (P 0 Drawer QQ) 
Santa Barbara, CA 93102 
Olcy Attn DASIAC 

General Research Corp. 

Santa Barbara Division 
P O Box 6770 
Santa Barbara, CA 93111 
Olcy Attn B Alexander 

Higgins, Auld Association 
2601 Wyoming Bi d NE 
Albuquerque, NM 87112 
Olcy Attn J Bratton 
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III Research Institute 
10 W 35th Street 
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Olcy Attn Documents Library 
Olcy Attn R Welch 
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Information Science, Inc. 

123 U Padre Street 
Santa Barbara, CA 93105 
Olcy Attn W Dudziak 
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Olcy Attn R Jones (Unclas only) 

Martin Marietta Corp. 
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Olcy Attn G Fotieo 
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5301 Bolsa Avenue 
Huntington Beach, CA 92647 
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Altadena, CA 91001 

Olcy Attn W Green 

Nathan M. Newmark Consult 
Eng Services 
B106A Civil Eng Bldg. 
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